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By large scale Monte Carlo simulations it is shown that the stable fixed point of the SO (5) 
theory is either bicritical or tetracritical depending on the efi'ective interaction between the antifer- 
romagnetism and superconductivity orders. There are no fluctuation-induced first-order transitions 
suggested by e expansions. Bicritical and tetracritical scaling functions are derived for the first time 
and critical exponents are evaluated with high accuracy. Suggestions on experiments are given. 
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High- Tc superconductivity (SC) in cuprates is achieved by carrier doping into the half-fiUing insulator when the 
antifcrromagnetism (AF) is suppressed. The SO (5) symmetry between AF and rf-wave SC has been proposed by 
Zhang which should be achieved at the bicritical point where the two critical {second- order) lines merge Q]. The 
variation of the chemical potential from the bicritical point breaks the SO (5) symmetry into the SO (3) for AF and 
the U(l) for SC. The question on the asymmetry between AF and SC caused by the Hubbard gap has been resolved 
by introducing the projected S0(5) concept Reasons that the S0(5) symmetry can be realized in reality have 
been discussed successfully by Arrigoni and Hanke Although the SO (5) theory gives a natural explanation on the 
phase diagrams observed experimentally, there is still no direct observation of the bicritical point in cuprates, partially 
because of the difficulty of chemical preparation of samples. Observations favoring the SO (5) scenario are reported 
in other systems: a bicritical point in a set of organic superconductors [Q, and more recently a tetracritical point in 
heavy- fermion superconductors mainly from the shape of phase diagrams. Since experimental data are piling in a 
wide range of materials, it becomes important to test the SO (5) theory in real systems in an unambiguous way, such 
as measuring the critical exponents To make this approach feasible in practice, and to resolve more fundamental 

problems discussed below, better theoretical understanding about the bicritical phenomenon of the SO (5) theory is 
necessary. 

Theoretical studies on bicritical phenomena go back to more than two decades ago § . According to RG e expansions 
a bicritical point becomes instable when the degrees of freedom n exceeds a critical value, ric ~ 4 in three dimensions 
(3D), and a so-called tetracritical point with four critical lines crossing takes over Q. Since rt = 5 in the S0(5) theory 
is only slightly above ric, the stable fixed point is nontrivial. As a matter of fact, basing on e expansions Murakami and 
Nagaosa argued 1^ that the bicritical point occurs only exceptionally; instead one should observe either a tetracritical 
point or a triple point where two first-order lines merge into a single first-order line (see Ref. ||l0[| for a more general 
discussion about possible phase diagrams). Since the S0(5) symmetry is achieved only at the bicritical point, a full 
understand on its stability is of fundamental importance for an assessment of the S0(5) theory. 

The RG is not necessarily correct in 3D when e expansions are truncated at low orders. A particularly relevant 
example is the normal to SC transition, which is fluctuation-driven first order in d = 4 — e, but is in XY-like second 
order in 3D. In order to derive the correct fixed points of the S0(5) theory by a strong-coupling approach, we have 
performed large-scale, classical Monte Carlo (MC) simulations. The main results are as follows: (1) The bicritical 
(tetracritical) point is stable for repulsive (attractive) interactions between the AF and SC orders. There are no 
fiuctuation-induced first-order transitions suggested by e expansions; (2) Bicritical and tetracritical scaling functions 
are derived for the first time; (3) The critical exponents and the ratios between coefficients of the two critical lines 
are evaluated with high accuracy. Since these properties are universal for the SO (5) theory, microscopic details are 
irrelevant. This opens the way to test unambiguously the S0(5) theory by experiments. 

The Hamiltonian in the present study is pj] ] 




(1) 
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on the simple cubic lattice ||T^. The coupling is limited to nearest neighbors. The vector s(t) of three(two) components 
is for the local AF(SC) order parameter, and sf + = 1 The symmetry-breaking g field is related to energy 

gaps of magnons and charged pairs, and the chemical potential of electron by g — (A^ — Ac + /i) /2 The w field 
governs the interaction between the AF and SC orders. Since terms of sf and can be mapped into the g and 

w terms up to a constant, and higher order terms are irrelevant to the critical phenomena, the above Hamiltonian is 
general enough for investigation on the multicritical phenomena of the SO (5) theory [p^ . In the above Hamiltonian, 
w >, = , < corresponds to >, = , < in Ref. respectively. 

In a typical simulation process, we generate a random configuration of superspins at a sufhciently high temperature, 
and then cool down the system gradually using the standard Metropolis algorithm of the MC technique. The g and w 
fields are fixed in this simulation process. Typically, we use 50,000 MC steps for equilibration and 100,000 MC steps 
for statistics on each site in the system at each temperature. The system size is = 40^ with periodic boundary 
conditions. We have also simulated systems of Li^ = 60'^ for several parameter sets and confirmed that the finite-size 
effects do not change the results. 

In order to establish the scaling analysis and to make accurate estimates of critical exponents, let us start with 
w = 0, taking advantage of the exact bicritical value gh = 0. The temperature dependence of the helicity modulus 
1^,^, which is proportional to the superfluid density and thus serves as the long-range order parameter of SC, is 
shown in Fig. 1(a) for positive g fields. For negative g fields we cannot observe the SC long-range order. The critical 
temperature and the critical exponent defined in T ^ (1 — T /Tc)"" are evaluated for g ~ 0.5J as v = 0.666 ± 0.004 
and Tc = 0.9773 ± 0.0004 J/fcs as shown in Fig. 1(a), where data points in OSJ/ks <T < OMJ/ks are used 0. 
This estimate of the critical exponent is in good agreement with 1^2 = 0.669 ± 0.002 by RG for the pure 3D XY (n=2) 
model noting v = {D — 2)v2 where = 3. Therefore, at least for g > 0.5 J there is a sufficiently large critical 
region where AF fluctuations are irrelevant to the SC critical phenomenon. 

The temperature dependence of the Neel order parameter is shown in Fig. 1(b) for negative g fields. We cannot 
observe the AF long-range order for any positive g field. The Neel point and the critical exponent defined in = 
W V< T^r^'i > ^ {I- T/TnY^ are evaluated as = 0.363 ± 0.001 and Tn = 0.9357 ± 0.0002 J/fcs for 5 -0.5 J as 
shown in Fig. 1(b), where data points in temperature region 0.7 J/ ks < T < 0.92J/kB are used. This estimate is in 
good agreement with the value = 0.3645 ± 0.0025 by RG for the pure 3D Heisenberg (n=3) model [^. Therefore, 
at least for g < —0.5 J there is a sufficiently large critical region where SC fluctuations are irrelevant to the AF critical 
phenomenon. 

As seen in Fig. 1, the critical onsets of the order parameters below the corresponding critical points become steeper 
in both SC and AF regimes as the absolute value of the g field is reduced. In order to analyze the crossover phenomena 
accurately, we extend the scaling theory for divergent quantities above the critical lines | p9[ to order parameters below 
the critical lines. We hypothesize that the helicity modulus for T < Tc{g) where g > gt — should depend on 
temperature and the g field in a scaled way 

T(r; g) ^ A{g - g,r/^ x /[(T/T^ - l)/{g - g,)^/'^] (2) 
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with the critical exponent for correlation length for n = 5 and the crossover exponent (j). The scaling function 
f{x) should have the following properties: (i) f{x) ~ (— a;)"^^, as a; ~oo (i.e. g gb for T < Tt,), which covers the 
critical phenomenon at g = gi,; (ii) f{x) ~ (i?2 — , as x B2 with B2 defined in B2 x (g — gt)^^'^ ~ Tc{g)/Ti, — 1 
(i.e. T — > Tc{g)), which accommodates the singularities at Tc{g)^s for g > gi,. From the scaling ansatz (^), one should 
have 

T(r,; g)/T(T,; g') = [{g - g,)/{g' ~ h)^/^ (3) 

Therefore, the bicritical point can be estimated as the temperature where the ratio on the left-hand side depends 
merely on the ratio (.g — 5b)/(.g' — ^b) but not on g fields themselves. Meanwhile, the combination of exponents vs/cf) can 
be determined. From the temperature dependence of the ratio T{T; g) /T{T; g') displayed in Fig. 2(a) for [gi,g[] = 
[0.05 J, 0.10 J] and [32,32] = [0.10 J, 0.20J] with gi/g[ = gaM = 0-5 (.96 = 0), we estimate Tt ~ 0.8458 ± 0.0005 J/Zcb 
and 1^5 /(j) ~ 0.525 ±0.002. Each point in Fig. 2(a) is obtained by 3 million MC steps on each superspin in the system. 

The crossover exponent can be evaluated from the slopes ki of T{T; gi)/T{T; g[) and K2 of T {T ; g2) / T {T ; g2) 
with respect to temperature at the bicritical point Tf,: 

(j) = ln[(32 - gb)/{gi - gb)]/ ln(Ki/K2). (4) 

Fitting the temperature dependence of T(T; g)/T(T; g') up to the second order of deviation from the bicritical 
temperature as shown by the solid curves in Fig. 2(a), one obtains ki ~ 7.7301 ± 0.0804 and K2 — 4.6907 ± 0.0309. 
Therefore, we have = 1.387 ± 0.030 with the error bar estimated from those of the slopes. Finally, we arrive at 
1^5 = 0.728 ±0.018. 

Our estimates of the critical exponents should be compared with the results by the e expansions of RG [^0|j2l[ 
— 0.723 (up to the third order of e) and — 1.313 (up to the second order of e). The agreement between the two 

estimates on 1/5 is extremely fine. The present estimate on is expected to be of higher precision than that of the e 

expansions. 

The scahng plot of Eq.(|) for all data in Fig. 1(a) is displayed in the main panel of Fig. 2. A magnification around 
zero helicity modulus is given in Fig. 2(b), where the solid curve corresponds to the critical exponent 1^2 = 0.666 and 
i32 — 1/4. It is clear that both properties (i) and (ii) are satisfied very well. 

It is worthy to point out that, using the scaling property discussed above, the bicritical point and the associated 
bicritical exponents can be inferred from informations away from the bicritical point. 

The same scaling analysis can be performed for the Neel order parameter for g < in Fig. 1(b). Since the Neel order 
parameter suffers more serious finite-size effects, the plot similar to Fig. 2(a) fails to give accurate estimate of critical 
exponents. Alternatively, we evaluate the order parameter 7715 = -^\/< > + < Si > g = gb ~ by MC 

simulation as shown in Fig. 3(a). The critical point and critical exponent are evaluated as = 0.400 ± 0.002 and 
Tb = 0.8559±0.0003 J/fcs> using data points in temperature region 0.7 J/ks < T < Q.SAJ/ks- The value of the critical 
exponent agrees very well with /3 = 0.402 by the e expansions to second order [pO|. A fine scaling plot is displayed in 
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the main panel of Fig. 3. From Fig. 3(b), one reads B3 ~ 1/6 with B3 defined in ,63 x (g^ — .9)^^"^ = Ti^{g)/Ti, — 1. 
The ratio 82/8^ = 3/2 is given very accurately by the inverse ratio between degrees of freedom (discussed in Ref. 
in a different context) and should be a universal constant for the SO (5) theory. 
We then turn to w > for which the fixed point is nontrivial as discussed previously. Without losing generality, let 
us set w = O.IJ. We observe the SC long-range order for g > 0.012 J and AF one for g < 0.010 J as shown in Figs. 4(a) 
and 5(a). At g = O.OllJ either the AF or the SC order is obtained depending on the initial random configuration and 
cooling process, same as at g = (75 = for w — 0. As clearly seen in Figs. 4(a) and 5(a), all the phase transitions upon 
temperature reduction are second order, and thus there arc no fluctuation- induced first-order transitions suggested by 
e expansions. No coexistence between AF and SC is observed in the g — T plane, and the SC and AF long-range order 
parameters reach their full values at zero temperature as in Figs. 4(a) and 5(a) (see Figs. 6(a) and (b) for the case of 
coexistence). Therefore, the fixed point for w — O.IJ is a bicritical point at gb ~ O.OllJ and Ti, — 0.8458 J/ks (equal 
to that for w = because of the same exchange coupling). Furthermore, all the data in Figs. 4(a) and 5(a) collapse 
on the bicritical scaling curves in the main panels of Figs. 4 and 5 with the same bicritical exponents obtained for 
w — 0, except for some deviations observed in the Neel order parameter for g — O.OIOJ coming from the finite-size 
effect. We observe the universal ratio B2/B3 ~ 3/2 in the present case too. From these observations, we conclude 
that the phase diagram for w > is governed by a bicritical point. 

For the negative cubic anisotropy w = —O.IJ, we observe a tetracritical point at gt = —O.OllJ (of the same absolute 
value of gi, for w = O.IJ) around which AF and SC coexist in a narrow region in the phase diagram displayed in Fig. 
6. From the scaling analysis shown in Fig. 7(a) with [51,. gi] = [0.039 J, 0.089 J] and [32,. 92] = [0.089 J, 0.189 J] with 
(51 - gt)/{g'i ~ gt) = (m - gt)/{g'2 - gt) = 0.5, we obtain Tt = 0.8460 ± 0.0005 J/ZcB (same with n for w > within 
the error bar), utl(t)t ^ 0.529 ± 0.002, and ki ~ 7.8900 ± 0.0935 and K2 ^ 4.7451 ± 0.0250. Adopting the same process 
developed for w = 0, we arrive at vt =i 0.721 ±0.018, ~ 1.363 ±0.030, and ft ~ 0.399 ±0.001, where ft is evaluated 
from data in Fig. 8(a). The critical and crossover exponents for the bicritical and tetracritical phenomena are too 
close to each other to be differed numerically within the errors of estimation. It should be difficult for experiments 
too. Our estimates on the tetracritical critical and crossover exponent are more accurate than the e-expansion results 
vt — 0.632 and (pt — 1.155 to the first order of e [||J|]. Deviations from the scaling curves in the main panels of Figs. 
7 and 8 occur when the AF and SC orders coexist. The universal ratio B2/B3 = 3/2 is observed once again. 

p — T phase diagram with the charge density p — (tf ) /2 can also be composed such as the one shown in Fig. 1(c) 
for w = 0. The two lines below the bicritical point are given by p±o{T) = ]iu\g^±{)p[T;g). The phase surrounded 
by these two curves corresponds io g — gh — Q- In this phase SC and AF coexist in a macroscopic level, while in 
a mesoscopic or microscopic level there can be phase separations, such as stripes. The discontinuous jump in the 
charge density at g = = and T < Tb is described numerically by p+o{T) — p-o{T) ~ (1 — T/Th)^" . We obtain 
ft = 0.80 ±0.08 using the scaling relations ft = 2 — Q!5 — ||] and — 2 — Dv^ — —0.184 ±0.054. As temperature is 
reduced at g = gb = 0, the charge compressibility diverges as dp/dg ^ (T — Ti,) with 7p ~ 20— 2 + 05 ~ 0.59±0.11. 
This divergent charge compressibility might be responsible for the peculiar transport properties observed below the 
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pseudogap temperature in cuprates. The shapes of the critical hnes in the p — T phase diagram Fig. 1(c) are also 
consistent with the behaviors ±(p — 0.2) ~ {Tc^n /Tb — 1)^'' as predicted by the scaling theory 

Finally, a suggestion for experimental investigations is made. As revealed in our study, by means of the scaling 
properties it is able to infer the bicritical (or tetracritical) point and the associated critical phenomena even in cases 
that sample preparations are difficult in the merging region of the AF and SC orders, such as high-Tc cuprates. 
Therefore, systematic measurements on the London penetration depth (which is related to the SC order parameter 
discussed in the present study by A ~ 1/ a/T) for a series of dopings near the boundary to the AF phase (see Ref. ]2^] 
for optimal doping) should be able to clarify the symmetry governing the high-Tc superconductivity. 
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Figure Captions 

Fig. 1: g — T phase diagram with a bicritical point for w = 0. Inset (a): temperature dependence of the helicity 
moduh; inset (b): that of the Neel order parameters; inset (c): p — T phase diagram. 

Fig. 2: Scahng plot for the hchcity moduh shown in Fig. 1(a) for w = with Tf, = 0.8458 J/fcij, = 0.728, = 1.387, 
and j/2 = 0.666. Inset (a): scaling plot based on Eq. (|^); inset (b): magnification of the main panel around zero 
helicity modulus. 

Fig. 3: Scaling plot for the Neel order parameter shown in Fig. 1(b) for w — Q with = 0.8458J/A;b, 13^ = 0.400, 
4> = 1.387, and (3^ — 0.363. Inset (a): temperature dependence of 5-component superspin order parameter; inset (b): 
magnification of the main panel around zero Neel order parameter. 

Fig. 4: Scaling plot for the helicity moduli shown in inset (a) for w ~ O.IJ with = O.OllJ and the same and 

critical exponents in Fig. 2. Inset (b): magnification of the main panel around zero helicity modulus. 

Fig. 5: Scaling plot for the Neel order parameter shown in inset (a) for w = O.IJ with gt, — O.OllJ and the same 

and critical exponents in Fig. 3. Inset (b): magnification of the main panel around zero Neel order parameter. 

Fig. 6: g — T phase diagram with a tetracritical point for w = —O.IJ. Inset (a): temperature dependence of the 

helicity moduli; inset (b): that of the Neel order parameters. 

Fig. 7: Scaling plot for the helicity moduli shown in Fig. 6(a) for w = —O.IJ with Tt — 0.846J/A:b, i^t = 0.721, 
(j)t = 1.363, and 1^2 = 0.666. Inset (a): scaling plot based on Eq. inset (b): magnification of the main panel 
around zero helicity modulus. 

Fig. 8: Scaling plot for the Neel order parameter shown in Fig. 6(b) for w = —O.IJ with Tt = O-SAGJ/ks, Pt = 0.399, 
(j)t = 1.363, and = 0.363. Inset (a): temperature dependence of 5-componcnt superspin order parameter; inset (b): 
magnification of the main panel around zero Neel order parameter. 
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